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System for Identifying Patient Response to Anesthesia Infusion 
Background of the Invention 

FIELD of the invention 

This invention relates generally to systems for monitoring patients during 
infusion of anesthesia, and more particularly, to a model structure that captures basic 
characteristics in patients' responses to anesthesia and surgical stimulation, such as BIS- 
based responses. 

description of the related art 

The goals of general anesthesia are to achieve hypnosis, analgesia, and patient 
immobility simultaneously throughout a surgical operation while maintaining the vital 
functions of the body. One of the most critical tasks of anesthesia is to attain an 
adequate anesthetic depth. 

Traditionally, anesthesia depth has been deduced indirectly from physiological 
signals, including autonomic responses such as heart rate, blood pressure, tearing, and 
patient's movements in response to surgical stimulations. The BIS monitor, which is 
based on the bi-spectrum level of an EEG signal, provides a viable direct measurement 
of anesthesia depth. The BIS index ranges from 0 to 100, where 0 corresponds to a flat 
line EEG, and 100 is a fully awakened state. Deep sedation is present at 60 and below, 
where the patient does not respond to verbal stimulus and has low probability of explicit 
recall. As a result of this technological advance, th^ have been many attempts to apply 
control methodologies to assist or automate the drag infusion. These preliminary studies 
on computer-aided drug infusion have been developed with the use of simple control 
strategies such as fuzzy logic or lookup tables. 

To facilitate further control strategy development for improved anesthesia 
performance, it is desirable to develop representative models of BIS responses to dmg 
infusion. Such models will allow more substantial and faster development and testing 
of new control, signal processing and decision methodologies. Also, they will provide 
a non-risky platform to test performance, robustness, and safety under extreme 
conditions and rare events. A satisfactory modeling method for this application must 
address several unique and challenging issues. Unlike electrical, mechanical, and 
chemical processes that are often repeatable and data rich, patients differ dramatically 
in metabolism and pre-existing medical conditions, as well as their responses to the 
contemplated surgical procedures. Consequently, individualized models must be 
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established for each patient with limited patient information and data points. The expert 
knowledge of anesthesiologists plays critical roles in predetermining patient response 
characteristics. Thus, there is a need for a model thait is sufficiently simple for easy 
clinical utility, suitable for development and testing of control methods for improved 
performance, and capable of incorporating expert knowledge when the data are 
insufficient. 

The research effort to develop computer-aided drug infusion systems has been 
both intensive and extensive since the early 1950s. The control methodologies 
employed include simple control {e.g., programmed control, relay control, and PID 
control), adaptive control self-tuning control, self-organizing control, and dual 
mode controller), and intuitive or intelligent control (e.g., fuzzy control, neural network 
control, and expert system-based control). To measure anesthesia depth, many indices 
from the EEG signals have been proposed, including the median frequency, spectral 
edge frequency, and auditory evoked potentials. More recently, there have been 
provided in the art BIS monitors having claimed medical and economic benefits 
including reduction in hypnotic drug usage, earlier awakening, faster time to meet 
post-anesthesia care unit discharge criteria, better global recovery score, and more 
accurate assessment of awareness risk. 

It is a great challenge to characterize mathematically a patient's response to drug 
infusion. As a result of large deviations in the aforestated physical conditioning, age, 
metabolism, pre-existing medical conditions, and responses to surgical procedures 
among various patients, there is a demonstrable high non-linearity and large variation 
in their responses to drug infusion. Physiology-oriented models are usually too 
complicated to establish individually using limited clinical data from a patient. On the 
other hand, anesthesiologists have been administering drug infusion successfully with 
only lintiited information on patients, such as weight and medical condition(s). The 
control strategies of an experienced anesthesiologist are based substantially intuitively 
on basic characteristics, such as the sensitivity of the patient to a drug infusion. There 
is, therefore, a need for a system that provides objective assistance to the 
anesthesiologist, thereby reducing the medical professional's reliance on subjective 
criteria. 

It is, therefore, an object of this invention to provide a system that provides 
predictive information of a patient's response to anesthesia and surgical stimulation. 

It is another object of this invention to provide a system that enables predictive 
control to compensate for surgical stimulation. 
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It is also an object of this invention to provide a system that provides visual 
indication of the predicted impact of a drug infusion decision on anesthesia depth and 
corresponding physiological variable. 

It is a further object of this invention to provide a system that provides indication 
5 to an anesthesiologist of impact to the status of a patient that is anticipated from surgical 

stimulation from predetermined aspects of a surgical operation, such as incision and 
closing. 

It is additionally an object of this invention to provide a system that signals to 
an anesthesiologist a warning of impending undesirable or critical patient conditions. 
10 It is yet a further object of this invention to provide a system that utilizes real- 

time data in combination with the anesthesiologist's assessment to identify 
predetermined characteristics that are particular to a patient. 

It is also another object of this invention to provide a system that facilitates 
optimization of drug dosage to achieve a desired patient status during surgery. 

15 Summary OF THE IPTVENTION 

Hie foregoing and other objects are achieved by this invention which provides 
a knowledge-based and control-oriented modeling approach. The models retain only the 
key characteristics of patient responses that are essential for control strategy 
development and can be determined by either data or expert knowledge. The existing 

20 knowledge of anesthesiologists is abstracted and analyzed from their manual control 

methodologies. Time delay, response speed, and drug sensitivity turn out to be the key 
properties in this regard. This understanding is then used to develop a basic but flexible 
model structure. Its validity and clinical utility are verified by actual clinical data. 

Although feedback control has been attempted in computer-aided or automated 

25 drug infusion, it will be shown herein that feedback control alone cannot avoid large 

fluctuations in BIS values when significant surgical stimulation is imposed. This 
drawback is largely due to the time delays in a patient's response to drug infusion. In 
this regard, predictive control becomes highly desirable. In predictive control, an early 
waming of predictable surgical stimulation such as incision will be sent to the control 

30 module. Based on an established model on the patient's response to the dmg infusion 

and to the surgical stimulation, an advance adjustment of drug infusion rates (feed 
forward control) can be administered to compensate for the impact of the stimulation. 
As a result, the BIS value fluctuations can be greafly attenuated beyond the capability 
of a feedback system. 
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In accordance with a &st method aspect of the invention, there is provided 
method of assisting a human expert in reducing predictable variations in the depth of 
anesthesia during the administration of a medical anesthesia drug to a patient. The 
method includes the step of solving: 

Xy X2 



where the coefficients Cj, C2, Q, as well as the time periods ^ (initial time delay after 
drug infusion) and Tp (time constant representing speed of response) are initiated by 
assessment of a human expert. 

In one embodiment of the invention, the human expert performs the step of 
10 assigning a relative value between 1 and 10 to represent the patient's response to 

infusion of the anesthesia drag, where 1 represents the slowest and 10 represents the 
fastest. 

Typical set points are selected to be approximately ^50, X2 » 100, andx^ ^ 150. 

In accordance with a further method aspect of the invention, there is provided a 
IS method of determining a model that corresponds to a predicted response of a patient to 

anesthesia drag delivery. The method includes the steps of: 

first determining an initial time delay ^ aft^ drag infusion for the patient; 

second determining a time constant representing speed of response of the 
patient; and 

20 third determining a nonlinear static function representing the sensitivity of the 

patient to a dosage of the anesthesia drag at steady state. 

In one embodiment of this further method aspect of the invention, the steps of 
first, second, and third determining are implemented in a Weiner structure. In a further 
embodiment, the steps of first, second, and third determining are implemented in a 

25 Hanmierstein structure. 

In accordance with a system aspect of the invention, there is provided a system 
for determining a predicted response of a patient to the administration of an anesthesia 
drag. The system includes a first memory for storing patient dynamics information 
relating to the infusion of a bolus dosage of anesthesia drag, the first memory having a 

30 first output for producing a first output signal corresponding to a first anesthesia level. 

There is additionally provided a second memory for storing patient dynamics 
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information relating to the infusion of a titrated dosage of anesthesia drug, the second 
memory having a second output for producing a second output signal corresponding to 
a second anesthesia level. A third memory stores patient dynamics information relating 
to the patient's predicted response to events of surgical stimulation, the third memory 
having a third output for producing a third output signal corresponding to an anesthesia 
effect level. A signal combiner arrangement receives the first and second output signals 
and the anesthesia effect level, and produces at an output thereof a combined anesthesia 
effect signal. A limiter is coupled to the output of the signal combiner for establishing 
maximum and minimum values of the combined anesthesia signal. Finally, a virtual 
anesthesia monitor produces an anesthesia value responsive to the combined anestiiesia 
signal* 

In one embodiment of this system aspect of the invention, the first, second, and 
third anesthesia levels correspond to respective BIS levels, the anesthesia effect level is 
a BIS level, and the combined anesthesia signal is a combined BIS level signal. In some 
embodiments, the virtual anesthesia monitoris a virtual BIS monitor for producing aBIS 
value responsive to the combined BIS signal. 

In a further embodiment of this system aspect of the invention, there is provided 
a source of known unpredictable disturbances for producing an unpredictable 
disturbances signal, and the signal combiner arrangement is arranged to receive the 
unpredictable disturbances signal and the combined anesthesia effect signal is responsive 
to the unpredictable disturbances signal. 

Brief DEscRiFnoN of the Drawing 

Comprehension of the invention is facilitated by reading the following detailed 
description, in conjunction with the annexed drawing, in which: 

Fig. 1 is a block and line representation of a patient model structure in 
accordance with the present invention; 

Figs. 2a and 2b are graphical representations that are useful in characterizing a 
titration model aspect of the invention; 

Fig, 3 is a block and line representation of a titration model stracture; 

Fig. 4 is a is a graphical representation that shows the relationship between a 
scaling factor and a BIS value; 

Figs. 5a to 5c are graphical representations of actual patient responses; 

Fig. 6 is a graphical representation of scaled surgical stimulation levels; 

Fig. 7 is a graphical representation of a drug sensitivity function (titration); 
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Fig. 8 is a graphical representation of a drug sensitivity function (bolus); 
Fig. 9 is a graphical representation of a stimulation sensitivity function; 
Figs. 10a to lOd are graphical representations of patient model responses; 
Fig, II is a block and line representation of a feedback only control arrangement; 
5 Fig. 12 is a block and line representation of a PID controller; 

Figs. 13a to 1 3e are graphical representations showing the quality of the feedback 
control; 

Fig. 14 is a block and line representation of a control arrangement; 

Fig. 15 is a block and line representation of a predictive and feedback control 
10 arrangement; 

Figs. 16a to 16e are graphical representations showing the control quality; 

Fig. 17 is a block and line representation of a system configuration showing 
certain ones of the modules; 

Fig. 18 is a block and line representation showing the core learning functions; 
IS Figs. 19a and 19b are graphical representations showing a patient' s BIS response 

to propofol titration; 

Figs. 20a to 20h are graphical illustrations that demonstrates the comparisons 
between recursive algorithms; 

Fig. 21 is a graphical illustration that shows the result of block recursion of a 
20 simulation; 

Figs. 22a to 22d are graphical representations showing a moving average model 
of order 40; 

Figs. 23a to 23d are graphical representations showing a moving average model 
of order 200; 

25 Figs. 24a to 24c are graphical representations showing an ARMA model of 

order I; 

Figs. 25a to 25e are graphical representations showing an ARMA model of 
order 40; 

Figs, 26a to 26e are graphical representations showing an ARMA model of 
30 order 10; 

Fig. 27 is a representation of typical assessment mappings; 

Fig. 28 is a representation of expert knowledge of drug impact on the outcome; 

Figs. 29a to 29c illustrate a two-step recursive estimation; 

Figs. 30a to 30d illustrate a patient model response; and 
35 Figs. 31a and 3 lb illustrate a real-time identified model. 
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Detailed Description of the Invention 

Fig. 1 is a block and line representation of a patient model structure in 
accordance with the present invention. The model contains three building blocks, each 
with a similar structure, as shown in Fig. 1. The first two blocks represent the patient 
5 response to propofol titration and bolus injection, respectively. The third block 

characterizes the patient response to predictable stimulation such as incision, closing, 
etc. For concreteness, it is assumed that propofol is the anesthesia drug although the 
model structure is generic and valid for other anesthesia drugs. 

1. MODEL STRUCTURES OF PATIENT DYNAMIC RESPONSES TO 
10 DRUG INFUSION AND STIMULATION 

The Ttfration Model 

Anesthesia drag titration, illustratively propofol titration, is administered in this 
specific illustrative embodiment of the invention by an infusion pump (Abbott/Shaw 
LifeCare™ Pump Model 4), The dynamic response between the titration conmriand to 
15 the pump and the drug infusion rate at the needle point is represented by a first-order 

dynamic delay: 

1 

7]^+ 1 

where Ti is the time constant of the infusion pump. T, typically ranges from 0.2 to 2 
seconds and can be predetermined from prior experiment data. 

The key component of the present model is the patient dynamics. Although the 

20 actual physiological and pathological features of the patient will require models of high 

complexity, for control purposes it is not only convenient but essential to use simple 
models as long as they are sufficiently rich to represent the most important properties 
of the patient's response. Understanding the information used by anesthesiologists in 
infusion control, the patient's response to propofol titration is characterized with only 

25 three basic components: (1) Initial time delay after drug infusion. During this time 

interval after a change of the infusion rate, the BIS value does not change due to time 
required for drags to reach the target tissues and to complete volume distribution. This 

is represented by a transfer function e"^''^ ; (2) Time constant representing the speed 

of response. This reflects how fast the BIS value will change once it starts to respond. 
30 This is represented by the transfer function: 
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1 

T,s+V 

and (3) A nonlinear static function representing sensitivity of the patient to a drug 
dosage at steady state. This is represented by a function or a look-up table/ 

Figs. 2a and 2b are graphical representations that are useful in characterizing a 
titration model aspect of the invention. The meanings of the three above-mentioned 
basic components are depicted in Figs. 2a and 2b. Consequently, a model structure for 
titration is shown in Fig. 3. 

Fig. 3 is a block and line representation of a titration model structure. The top 
portion of the figure is the Wiener structure in which the nonlinear function follows the 
linear part; the bottom portion of the figure is the Hammerstein structure in which the 
nonlinear function precedes the linear part. Both can used in this application. For 
concreteness, we will use the Wiener structure in this paper. The delay time T^, time 
constant Tp and the sensitivity function of the model vary dramatically among patients 
and must be established individually for each patient. Model development using clinical 
patient data will be detailed hereinbelow. 

Bolus Miection Model 

Bolus injection inputs a large amount of drag in a short time interval. The bolus 
model has a structure similar to the titration model. However, due to the highly 
nonlinear nature of the patient response, it has been observed that the patient response 
to bolus injection demonstrates diifferent model parameters. The dynamic part is 
represented by a transfer function 

Data analysis indicates that bolus injection response has a slightly shorter delay time due 
to an increased speed for drug to reach the target tissues. However, its time constant is 
much larger. Also, its residue impact on the BIS values is small at steady state. This 
results in a very different sensitivity function, which will be detailed in the next section. 
Surgical Stimulation 

One of the most important control tasks in infusion control is to predict the 
impact of surgical stimulation on the BIS levels. Surgical stimulation is a main cause 
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of large fluctuations in anesthesia depth during a surgical procedure. Early prediction 
of the impact will allow advance drug adjustment before a feedback signal from the BIS 
monitor activates control actions. Typical surgical stimulation includes initial 
preparation, incision, closing, etc. Such stimulations cause large fluctuations in BIS 

5 values and must be compensated by adjusting anesthesia drugs. An experienced 

anesthesiologist understands the level of such stimulation and adjusts drugs in advance 
to minimize its impact. 

Surgical stimulation has an opposite impact on the BIS values to the drug 
infusion. The higher the stimulation, the higher the BIS value. The sensitivity of the 

10 BIS values to stimulation decreases with deepening of anesthesia, Le. , lower BIS values. 

Unlike drug infusion whose rates can be controlled and measured, stimulation is not 
controlled by anesthesiologists and cannot be accurately measured. On the other hand, 
unlike unpredictable disturbances or noises, significant surgical stimulation is 
predictable in advance and can be roughly characterized by their severity in affecting 

IS BIS values. 

Classification levels from 1 to 5 to represent the severity of the stimulation: 1 
being a rmnov stimulation and 5 the most significant stimulation. For control purposes, 
these rough characterizations are clinically feasible and sufficient for predictive control, 
which will be described in detail hereinbelow (see. Predictive Control). Patient 

20 responses to surgical stimulation are represented by the third model block. The block 

contains also a time delay Tg and a transfer function: 

1 

to reflect its dynamics and a lookup table for its impact on the BIS levels at the steady 
state. Since a patient's BIS sensitivity to stimulation depends on hypnosis state, a scaling 
factor a is introduced, a takes values between 0 and 1. The surgical stimulation is first 
25 scaled by: 

Stimulation Level x a = Scaled Stimulation Level 

The scaled stimulation is then used to predict the impact on the BIS measurements. A 
typical curve of a as a function of the BIS value is depicted in Fig. 4. 

Fig. 4 is a is a graphical representation that shows the relationship between a 
30 scaling factor and a BIS value. Information on an upcoming stimulation will be used as 
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a pre-waming to provide a predictive control so that the impact of the stimulation can 
promptly be compensated by early tuning of drug infusion. 

BIS MONITOR DYNAMICS 

The BIS monitor measures the EEG signals and performs bi-spectrum statistical 
5 analysis. Due to the processing time and data windows required in this procedure, the 

BIS monitor introduces a pure time delay (from processing time) and a dynamics delay 
(related to the size of data windows). Its relevant dynamics for effecting control can be 
represented by a system transfer function: 

^- V — \ — 

cascaded with a nonlinear linniter with lower limit 0 and upper limit 100 (defining the 
10 BIS range). BIS monitor dynamics can be established off-line, independent of surgical 

procedures. 

Other Disturbances 

There are many sources of disturbances that will also affect the BIS readings. 
These include unpredictable stimulation to the patient, electrical noises, sensor 
15 attachment movements, environment noises, etc. These disturbances are unpredictable 

and cannot be compensated by early predictive control. Their impact on BIS levels will 
be minimized by a feedback action. 

MODEL VERIFICATION AND DEVELOPMENT BY CLINICAL DATA 
Data Collection Protocol 

20 To verify the utility of the model structure and develop individual patient 

models, clinical data are collected, after an approval from the HUB board. In the clinical 
data collection process, standard medical care is given to the patients without any 
deviations. A typical protocol is described in the following: The BIS monitor is 
applied. After the monitor sustains a base level of the BIS values, versed (a drug) and 

25 fentanyl are given for inducing preoperative sedation and controlling a BIS level higher 

than 90. After the patient is brought into the operating room, other physiological 
parameters will be monitored, such as blood pressure, EKG, pulse oximeter, and 
end-tidal CO2 (after administering oxygen). Sedation will start with 1 mg versed and 
bolus propofol injection of 20 mg IV, followed by propofol infusion of initial rate 25 

30 mg/kg - min to control the BIS level to the range of 70 - 85. Before surgical incision, 50 

mg of fentanyl is given and propofol infusion rate is adjusted to control the BIS level to 
50 - 65, After incision, fentanyl and propofol are continuously adjusted to maintain the 
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BIS level at 60 - 75 until the end of the surgery. When the surgical closing starts, the 
BIS level is adjusted to 70 - 85 until the end of anesthesia drug administration. 

Actual drug infusion rates and BIS trajectories vary with surgical procedures. 
Drug infusion rates are manually recorded. BIS trajectories are recorded in the monitor 
5 and transferred to a laptop computer using the standard HypeiTerminal software. The 

data are then translated into Matlab*™ data files for further data processing. 

Clinical Data and Analysis 

One typical data set is used in the following subsections to illustrate the 
modeling process and verification results. In this case, the anesthesia process lasted 

10 about 76 minutes, starting from the initial drug administration and continuing until last 

dose of adnoinistration. Propofol was used in both titration and bolus. Fentanyl was 
injected in small bolus amount three times, two at the initial surgical preparation and one 
near incision. Analysis shows that the impact of fentanyl on the BIS values is minimal. 
As a result, it is treated as a disturbance and not explicitly modeled in the present 

15 example. The drug infusion was controlled manually by an experienced 

anesthesiologist. The trajectories of titration (in fdg/sec) and bolus injection (converted 
to g/sec) during the entire surgical procedure were recorded, which are shown together 
with the corresponding BIS values in Figs. 5a to 5c, which are graphical representations 
of actual patient responses. 

20 The patient was given bolus injection twice to induce anesthesia, first at r = 3 

minute with 20 mg and then at / = 5 minute with 20 mg. They are shown in the figure 
as 10000 iJLg/sec for two seconds, to be consistent with the titration units. The surgical 
procedures were manually recorded. Three major types of stimulation were identified: 
(1) During the initial drog administration (the first 6 minutes), due to set-up stimulation 

25 and patient nervousness. This is characterized as a level 3-3.5 stimulation by the 

anesthesiologist. (2) Incision at r = 45 minutes for about 5 minutes duration. This is 
characterized as a level 5 stimulation. (3) Closing near the end of the surgery at r = 60 
minutes. This is characterized as a level 3.5 stimulation. Due to different BIS values 
in these time intervals, the scaling factor a has different values under these three types 

30 of stimulation. The scaled stimulation level trajectory is shown in Fig. 6, which is a 

graphical representation of scaled surgical stimulation levels. 
Model Development 

The data from the first 30 minutes are used to determine model parameters and 
function forms. The infusion pump model and the BIS monitor model can be derived 
35 or identified with limited impact from real data. For the systems used in this study, the 
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time constant for the infusion pump model is estimated as = 0.5 second. The delay 
time and time constant of the BIS monitor depend on the setting. For a short window 
selection (for fast response), these parameters are estimated as = 10 seconds and 
Tf, = 10 seconds. 

For estimating the parameters in the titration block, the data in the interval where 
the bolus and stimulation impact is minimal (between / = 10 to r = 30 minutes) are used. 
By optimized data fitting (least-squares), tiie estimated parameter values ^ = 10 second 
and Tp = 90 second are derived. The data-generated sensitivity function is shown in Fig. 
7 which is a graphical representation of a dmg sensitivity function (titration). 

The first interval of data is then used to estimate the bolus model and stimulation 
model. The data are first processed to remove the impact of titration by using the 
titration model developed above. The processed data are then used to determine model 
parameters and function forms. The resulting parameters are (all in seconds): ttoto,= 10, 
r^to* = 250, r, = 12, r, = 3, X, = 4. The sensitivity functions are shown in Fig. 8, which 
is a graphical representation of a drug sensitivity function (bolus), and in Fig. 9, which 
is a graphical representation of a stimulation sensitivity function. 

Verification 

The actual BIS response is then compared to the model response over the entire 
surgical procedure. Comparison results are illustrated in Figs. 10a to lOd, which are 
graphical representations of patient model responses. Here, the inputs of titration and 
bolus are the recorded real-time data. The model output represents the patient response 
very well. In particular, the model catches the key trends and magnitudes of the BIS 
variations in the surgical procedure. This indicates that the model structure contains 
sufficient freedom in representing the main features of the patient response. Also, the 
impact of su]:g;ical stimulation is captured by the model. These features are of essential 
importance in designing control strategies for computer-aided drag infusion. 

FEEDBACK CONTROL 

Based on the patient model established hereinabove (Model Stmctures of Patient 
Dynamic Responses to Drug Infusion and Stimulation), a control strategy is described 
for maintaining BIS levels and to provide robustness. The control consists of two 
elements: (1) a feedback part for fine tuning titration to compensate BIS deviations due 
to drifting of patient dynamics, unpredictable disturbances and noises, model errors due 
to simplification and non-linearity; and (2) a predictive control to provide fast 
compensation for large but predictable stimulation. 
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To justify the utility of this control structure, there is first presented a discussion 
on feedback control and its performance. 
Feedback Control 

Feedback control relies on the measured BIS signal to activate control actions 
in order to attenuate BIS deviations due to disturbances and stimulation. Jn feedback 
control, even the predictable stimulation will be considered as disturbances to the 
system. The structure of feedback-only controllers is depicted in Fig- 11, which is a 
block and line representation of a feedback only control arrangement. 

It has been learned from inultiple studies of feedback controllers for use in 
anesthesia drug control that the highly nonlinear nature of the patient response precludes 
the direct application of linear control theory. Some controller types that can be tuned 
for nonlinear systems include PBD, fuzzy, etc. In the present embodiment, a typical PID 
controller employed to investigate feedback performance. 

The PID controller is illustrated as a block and line representation in Fig. 12. 
The parameters Kp, and are design parameters. It is well understood that choosing 
these parameters is subject to fundamental tradeoffs. For instance, large Kp will result 
in oscillations, while small Kp will result in a large control error. Similarly, large will 
amplify noises. Fine tuning of the parameters along these principles has resulted in the 
control parameters indicated in the figure. 

Performance Analysis and Fundamental Hmtfations 

Figs. 13a to 13e, which are graphical representations show the quality of the 
feedback control performance. A careful examination of the BIS trajectory reveals that 
the feedback control is very effective in attenuating BIS deviations due to random noises 
(see the portion of Fig. 13d where no large surgical stimulation exist). This observation 
indicates that feedback control design for attenuation of random or small disturbances 
in this application is not particularly difficult, even though the system is highly 
nonlinear. This is consistent with other reports on BIS-based closed loop anesthesia 
control. However, in the presence of large surgical disturbances, performance 
deteriorates very rapidly. This is clearly reflected at the time intervals where incision 
and closing occur. This phenomenon is of generic nature. One explanation for this 
phenomenon is that the BIS values respond to surgical stimulation faster than to drug 
infusion. As a result, there is a feedback loop time in which the influence of the 
stimulation is not compensated, regardless the control algorithms. Furthermore, the 
conmionly used method of relying on the BIS output speed to accelerate control actions 
in the event of stimulation cannot be effectively used in this application. This is due to 
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the fact that the BIS measurements are derived from EEG signals, which contain a lot 
of high frequency noises. Acceleration of feedback will cause amplification of noises. 
Consequently, improved control strategies are needed beyond pure feedback control, as 
will be described below. 

5 PREDICTIVE CONTROL 

Design and Performance 

It is observed from clinical data that it typically takes a delay time of 5-20 
seconds before the BIS value changes after a modification to the titration rate. The BIS 
value changes gradually thereafter with a time constant 80-200 seconds. On the other 

10 hand, a surgical stimulation impacts the BIS levels much faster. Consequently, a 

surgical stimulation will inevitably cause large deviations on the BIS values if the 
control action depends solely on the BIS monitor. This will significantly compromise 
the quality and safety of the anesthesia care of the patient. Such large deviations are 
evident from the manual controlled BIS trajectory in Fig. 10. For example, the BIS level 

15 rises to 90 after incision, which is highly undesirable. From a control viewpoint the only 

remedy is to provide preventative drag infusion in prediction of a significant stimulation. 
This becomes the main goal of the predictive control block. 

The control structure is depicted at a high level in the block diagram of Fig. 14, 
which is a block and line representation of a control arrangement. It is noted that 

20 predictable stimulation is distinguished from noises. In accordance with the present 

control strategy, the adverse effects of noises are dinndnished by feedback control and 
predictable stimulation and axe compensated by predictive control algorithms. 
Information on predictable stimulation is sent in advance to the control module. 

Fig. 15 is a block and line representation of a predictive and feedback control 

25 arrangement showing a detailed control stmcture. For the patient characteristics under 

consideration, a simple feedback stmcture appears sufficient. Here, there is employed 
a basic proportional feedback (Kp = 10, the same as in the feedback-only control) 
together with a filter that compensates infusion pump dynanwcs. The predictive control 
block is significantly more involved. It is designed on the basis of optimal stable 

30 inversion of the patient dynamic models. In principle, the desired drug infusion rates 

under a given stimulation can be estimated by inverting the patient model. However, 
such an inversion will result in a non-causal and unstable system, and hence is not 
acceptable. The task here is to find a causal stable (nonlinear) system that can optimally 
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approximate the model inverse in the frequency band of interest. A computer program 
to find such an optimal system is then devised. 

To understand the potential benefits and limitations of the control strategy, we 
applied it to the patient model. The BIS trajectories from manual control are compared 
to those from feedback and predictive control. The simulation results are presented in 
Figs. 16a to I6e, which are graphical representations showing the control quality. Due 
to its predictive capability and robustness in compensating BIS deviations continuously, 
BIS levels can be maintained much more steadily under the computer control strategy. 

It should be emphasized, however, that if stimulation is not well predicted, then 
control quality will deteriorate. On the other hand, even incomplete information on the 
upconoing stimulation can be used to improve BIS fluctuations. 

Robustness against Modeling Errors 

Since a patient model characterizes only the key features of the patient response, 
it is only an approximation of the actual patient response. As a result, it is essential that 
the control strategies remain functional under modeling errors and inaccuracy in defining 
stimulation levels and scaling factors. The inventors herein have tested the robustness 
of this control stmcture by introducing many variations of potential modeling errors and 
perturbed stimulation characterizations. The predictive plus feedback control structure 
provides satisfactory performance under these perturbations. The actual testing results 
are omitted here for simplicity. 

Modules and Learning Functions 

Fig. 17 is a block and line representation of a system configuration showing 
certain ones of the modules. As shown in this figure, a plurality of signals, that may, in 
certain embodiments of the invention, include the BIS monitor signal, a blood pressure 
signal, a heart rate signal, and other signals is delivered simultaneously to a Real-time 
Learning and Identification Module after signal synchronization and to a Safety 
Assessment Module that compares the signal values to predetermined safety 
specifications stored in a database. The result of the safety assessment is displayed on 
a Display Interface. As shown, the Display Interface includes graphical and numerical 
displays, as well as warning indicators and messages. 

The Real-time Learning and Identification Module additionally receives 
information responsive to the assessment by an expert {Le. , an anesthesiologist) that has 
been processed by a Knowledge Mapping Module, The output of the Real-time 
Learning and Identification Module is conducted to an Anesthesia Knowledge-based 
Wiener Model Structure, that functions as hereinabove described. 
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Data corresponding to a desired target status for the patient (not shown) is 
delivered to an Optimal Drug Dosage Calculation Module that also receives the data 
from the Anesthesia Knowledge-based Wiener Model Structure. The status of the 
Optimal Drug Dosage Calculation Module is displayed on the Display Interface, and the 
corresponding data is delivered to a Stimulation Compensation Module that also 
receives the data from the Anesthesia Knowledge-based Wiener Model Stmcture. The 
data from the Stimulation Compensation Module is displayed. In addition , the data from 
the Anesthesia Knowledge-based Wiener Model Stmcture is delivered to a Dmg Impact 
Prediction Module, the output of which is also displayed. 

Fig. 1 8 is a block and line representation showing the core learning functions and 
their interrelationship in a specific illustrative embodiment of the invention. As is seen 
in this figure, two forms of input data (z.e., expert assessment and BIS and physiological 
data) are processed and delivered to an iterative process. 

The prior expert assessment data is organized in accordance with a 
predetermined knowledge mapping relationship and combined with the dynamic and 
sensitivity processing of the Weiner model to form a body of prior information. This 
prior information is then delivered to the first stage of the iterative process where the 
sensitivity portion of the model map data is captured. Simultaneously, the BIS and 
physiological data is processed, scaled, and synchronized, and is also delivered to the 
first stage of the iterative process. 

In this specific illustrative embodiment of the invention, the output of the first 
stage of the iterative process is delivered to a second stage, where the parameters of the 
dynandc portion of the model map is updated by stochastic approximation and projected 
least squares processes. The resulting data then is delivered to the third stage of the 
iterative portion of the process where the dynanric portion of the model map is captured 
and an interim item of drug infusion data is developed The interim item of dmg 
infusion data then is delivered to the fourth stage of the iterative portion of the process 
where the parameters of the sensitivity portion of the model map are updated by 
stochastic approximation and projected least squares processes. This process then is 
repeated by real-time iteration where t = t+1, until an adequate Real-Time Individualized 
Patient Model is formed. 

The mathematical underpinnings of these processes is described below. 
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Expert Knowledge and Fast Identification 

The patient model must be individually developed to accommodate differences in 
patient characteristics, prcr-existing medical conditions, surgical types, and other de- 
viations. Since anesthesia infusion decisions must be made before substantial data 
are available, it is imperative that prior expert knowledge of an anesthesiologist is 
incorporated into the learning process. In this section, a basic identification method 
will first be developed. Some essential properties of the method will be established. 
It will then be applied to real dinical data in the subsequent sections. 

Systems and Expert Information 

The system under consideration has a Wiener structure with input u and output y: 

where G is a dynamic linear system^ f is a memoryless nonlinear function, and w 
and V are random noises. In our applications here, G and / represent the dynamic 
(or transient) and steady-state behavior of a patient's response to a drug infusion 
input, respectively. Simple examples of G include G{s) = {9 = T) or e""^*^q^ 
(0 = [r, T]'). Here and henceforth, for a vector or matrix denotes its transpose. In 
our algorith, for computerized data processing, G will be discretized and represented 
by an ARMA (Auto-Regression and Moving Average) model 

Xu = -airCjb.i OnXfe^n + h^Uk + biUjfc-i + - • + &n^ib-n + ^^A: = ^^'fc^ + t^fc 

where 

^ = [Oi, . . . , On, 6o> • • • > M- 
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A special case of the ARM A model above is an MA (moving average system), whtere 
ai, • . . , On and a;(* — 1), , . . , a;(4 - ri) disappear in the parameter and regressor, respec- 
tively.* 

The prior information on 9 will be a bounded set Q,u which is provided by prior 
5 expert knowledge. For example, suppose G{s) = yq^. The model parameter T 
has a clear physical meaning of response speed and allows direct inputs from expert 
knowledge. For instance, a patient with fast drug response will have a small T. A 
linguistic assessment 'fast response" can be translated by a pre-designed assessment 
Q mapping to a bounded set [T, T). This information will then be mapped to a boimded 
set on 9 after discretization. Typical assessment mappings look like the intervals in 
Fig. 27. This set will serve as prior information for the learning algorithm to be 
described in the next subsection. 



The nonlinear function / can also be parameterized. Usually, polynomials are 
used. However, polynomials of the form x^^^ + Cm-2a?'"""^ H h cia; + co is not conve- 
nient for incorporating expert knowledge. In many medical applications, / represents 
steady-state sensitivity: The impact of the input (drug amount) on the outcome. A 
physician's knowledge is often expressed by an approximate relationship between the 
input and outcome: Given typical input amounts a;i, . . . , the approximate ranges 
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of the predicted outcomes . . . i [y^j^ml- accommodate this knowledge, we 

replace . • . ,a;"*""^ by a collection of appropriate functions ^i{x) and introduce 
the following parameterization for /, Suppose that the prior information on / is 
given by m distinct mappings: rr^ (y^, yj, i = 1, . . . , m as shown in Figure 2. We 
5 parameterize / by 

m 

y = m = 'E^iUx) (1) 

i=l 

where ${(a;) is the mth order polynomial 

^i(x) satisfies ^i{xj) = 0, j i] ^i{xi) = 1. In other words, ^i{x) are indicator 
10 fxmctions on the set points . . .^Xm}- The unknown parameters will be denoted 
by © = [Cu - • . , CmY' Consequently, the prior expert information becomes precisely 
• the prior interval information on the parameters: Ci € [y^> yj, i = 1, . . • , m. 



The expert knowledge will serve as the prior information on the unknown param- 
eters 9 and 0. It will be generically expressed as 9 e^i and 0 G Here, Qi and 
^2 are bounded sets. In many applications, they are simply polyhedres. (See, Fig 28) 
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Basic Setup for Real-Time Learning Algorithms 

We propose a class of estimation algorithms, which are recursive and can be thought 
of as a relaxation procedure that uses most recently acquired information. The algo- 
rithms are a new twist of the least squares identification algorithms used in a wide 
variety of applications. Its essence is to use real-time data on the input {uk} and 
output {vk} as soon as they become available. One of the distinct features of our 
algorithms is: The identification task is accomplished in two steps. The first step 
requires mapping the sequenece of outputs back to get a sequence of "intermediate" 
output, which is used for estimating 9. Then the second step takes care of the estima- 
tion of 9. These estimates are obtained recursively and hence are suitable for on-line 
computing. The expert prior knowledge is used as constraints in the algorithms. The 
problem can be stated as follows: Let 9k and 0^ denote the estimates of the model 
parameters at the A;th step.^ We proceed to design an identification algorithm for 9k 
and Ga:, based on the prior information fii, 0,2, and past and present- observations 
on u and y, such that a least squares type cost function is minimized. The following 
issues are of importance in designing the identification algorithms. 

1. Real-Time Computational Bfiiciency: During real-time implementation, the 
data size N becomes larger and larger. Consequently, the computational bur- 
den of finding estimates can potentially become overwhelming when N is large. 
As a result, recursive algorithms^ which employ newly acquired data Uk and 
yis to update the estimate from 9k and to 9k-^i and 6jb+i, are much more 
preferrable than their off-line counter part. 

2. Nonlinearity: Due to the nonlinear model structure, the nmnerical solution may 
be difficult to obtain. Hence, it is desirable to find an approximate optimization 
problem that is convex or linear. 

3. Convergence: We need to establish convergence properties of the designed iden- 
tification algorithms. 

We update the estimates of 9 and 0 alternately and use the most recently infor- 
mation as soon as they become available. To do, we construct iterative procedures so 
that the estimates of 9 and 0 are computed iteratively one followed by the other. To 
be more specific, we introduce the two-step iterative identification setup below. The 
more specific algorithms will be given later. 

Two-step Procedure. Given the bounded sets J^i and which are given by use 
of thejexpert Jcnowledge, and performance indices ei{0) and 62(0), we start at A; = 0 
with ^0 and 0o being interior of and ^2, respectively. 

^The actual time interval between the kth and the {k + l)st steps depends on the sampling 
interval of the system. 
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Step 1: At time A;, suppose that the estimates 9k Skud Ok have been computed. To 
construct the next estiamtes, the essence of the least squares approach normally 
reques a cost function involves the the difference of observations and "weighted" 
regression be minimized. However, the output, or measurement, or observation 
5 Xfe is not available. One remedy is to map the available G^. back and obtain 

the "intermediate" output using the estimated nonlinear^sensitivity function 
with parameter 0^. Then use Uk and to obtain update ^a;+i- 

The update algorithm is a projected identification recursive scheme 

1. 9* = argmin^ei(5) subject to the constraint 0* € f^i. Denote the estimate 
jQ at the kth step without constraint by 9k and the constraint estimate by 

9k' At time + 1, 

2. If 9k+i e fii, then 9k^i ^ 9k^i. 

3. If 9k+i ^ Qu then 9k^i is the projection of 9k+i onto f2i, namely the closest 
25 element of Qi to Sk+i- 

Step 2: Now since 9k^i has obtained, we use Uk and Vk to update 0^+1 from ©a:. This 
is accomplished by minimize another cost function 62(6) again in the mean 
sqaures sense. 

The update algorithm is also of the projected least-squares identification: 

20 1. 0* = arg mine 62(0) subject to the constraint 0* € ^2. Denote the un- 

constrained and the constrained estimates at the (k + l)st step by 0jb+i 
and 02A:+i> respectively. 

2. If ^k^i G Q2, then ©jt+i = Ok+i^ 

25 3. If ^ ^2, then ©fc+i is the projection of Qk+i onto 1^2, namely the 

closest element of Q2 to ^ib+i. 

Recursive Algorithms 
Least Squares Algorithms 

Following the discussion in the previous section, we develop the two-step procedure 
30 here. Note that the statistics of {xk, (f>k) and (y*, may not be known, but their 
measurements are available. 

To begin, we first construct the unconstrained algorithms and obtain its con- 
strained counter part for 9k and then write the corresponding part for 0^. In what 
follows, to signifty the dependence of the cost function on k, we write it as ei^ki9)' 
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Consider the problem of minimi^g 

where Xk are obtained by a back mapping technique via 0^. Suppose that the matrix 

1=0 

is invertible with probability one for k large enough. Then, the minimizing 6 is given 
by 

1=0 

Equation (3) can be put into a recursive form by expanding 



0M = nil 

0 to get see [9, pp. 18-20]) 



.<=0 J 



= + ^];i^(i>k - <f>k^ • (4) 

Since a matrix is involved, the computation can be time consuming. To avoid the 
onerous task, by using the matrfac inversion lemma [9, pp, 18-20], we can rewrite the 
procedure recursively as 

15 ^k^i-^k i+j;^ii!'',^4>.\ _ (5) 



These are known as the recursive least squares formulas. 

To further simplify the procedure, taking a first-order (in ^j^^) expansion in (4) 
and (5) yields a linearized least squares approximation 

20 Ok^i = ek + ^k^(t>kdk + «'fc Vibl^fc - (f>k^kl /fiN 

nli = [i'n'M'k]n'^ 

The convergence of the sequences in (5) is assured by any conditions that assure the 
convergence of the least squares estimator. 
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Similaxly, using the cost function 

C2.fc(e) = i|;(yi-0'*i)% (7) 



and assuming that 



<=sO 



" {=0 



is invertible, we arrive at the recursive least squares procedure for estimating 6* 

l + ^H,-^*.' , (8) 

and the linearized least squares approximation 



Adaptive-Filtering-Type Procedure 

Replacing the matrix gain by a scalar sequence of step sizes leads to further 
simplification of (6). In fact, we obtain the following recursive algorithm 

' dk+i = % + enM^k - (l>k^kl (10) 
where ek is a sequence of real numbers satisfying 

Sk>Of €k-^0 as fc-^oo, and ^eAs = oo. (11) 

k 

Algorithm (10) is a stochastic approximation type procedure, and belongs to the 
family of adaptive filtering algorithms; see [6] for an updated accoimt on stochastic 
approximation algorithms. Commonly used step size sequences include = 0(1/A;), 
Ek = 0{l/kP^) for some 0 < 7 < 1, among others. Likewise, for 0, we construct the 
algorithm as 

ejfc+i = Gfc + Sk^klvk - *;0fc], (12) 
where we choose the same sequence {ek} as in (11). 
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Constrained Algorithms 

For the constrained algorithms, we use projection operators Hi and 112 to denote the 
projections onto fii and J^2» respectively. Recall that Cti = {0 : 6^ < 6^ < 9^ i = 
l,...,2n + l} and S72 0* < ©i, 1= respectively. Thatjs, we 

let the ith component of the estimates 0 and Q be confined to the interval [g^, 0i] and 
[G^, ©i], respectively. Similar to the previous cases, the constrained algorithms can 
be written as r r n 



Variants 

Comparing the least squares algorithms with that of the adaptive filtering algorithms, 
their differences are seen in the step size sequences involved. Consequently, the adap- 
tive filtering type algorithms have the advantage over that of the least squares al- 
gorithm since one need not warry about matrix manipulations for the step size se- 
quences. As far as the computational complexity is concerned, the real-valued step 
sizes are more preferrable. Nevertheless, it is also clear that the least squares algo- 
rithms use more information compared to the adaptive filtering procedures owing to 
the use of the matrix gain sequence. Therefore, it is expected that the least squares 
algorithm provided better performance owing to the use of matrix gain. One question 
is almost immediate. Can we construct an algorithm that has the advantage as that 
of the adaptive filtering procedures with reduced complexity, in the mean time, still 
provides good performance? In other words, we wish to construct algorithms that 
have advcmtages of both adaptive filtering type and that of the basic least squares 
procedure. The ideas of averaging becomes a favorable alternative, which was initially 
considered in [12, 14] and subsequently treated in [5, 20] for much more general noise 
processes; see also [1, 15, 20] for averaging applied to both iterates and observations. 
In what follows, we present two of such averaging procedures. Both of them give rise 
to asymptotically efficient schemes. We will show that the averaging algorithms so 
constructed yield asmptotic optimality. By asymptotic optimality, we mean that not 
only is the scaling factor the best one, but also the limit covariance is the smallest 
possible. 

Suppose that 0k — > 0* the true parameter with probability one. For the rate of 
convergence, we examine the sequence k^{6k — ^*)- We try to find such real num- 
ber a that the above rescaled sequence converges to a nontrivial limit in the sense 
of convergence in distribution. The scaling factor A;^ together with the asymptotic 
covariance then gives us the desired convergence rate. 
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Constant-stei>*size Algorithms 

In lieu of decreasing step-size sequences used in (13), we may consider algorithms with 
constant-step size. In fact, in various signal processing applications, one normally uses 
a constant step size rather than a sequence of decreasing step sizes. The rationale is 
5 that the use of a constant gain normally enables one to track slight variation of the 
parameter. In reality, the parameters involved in patient responses to anesthetic drug 
infusion are not fixed constants, but rather time varying. Thus, the tracking ability 
is important. For our problem, suggest the following constant-step algorithms with 
constraints. ^ _ 



Note that in (14), 6k and Gfc+i in fact depend on e and should have been written 
as 6% and O^+i, respectively. For notational simplicity, we have suppressed the e- 
dependence. To be able to get meaningful asymptotic results for algorithm (14), we 
need to let e — > 0. In the actual computation, one uses a small fixed constant e. 

Iterate Averaging 

The idea behind this scheme is: First one obtains a rough estimate via the use of 
larger than 0{l/k) step sizes. Then one takes an averaging of the iterates. One of 
the particular attractive features of the algorithm is that the use of large step size 
20 creats oscilations and force the iterates moving towards the target faster than the use 
of smaller step sizes. In addition, although the averaging is essentially ofi'-line, it can 
be executed recursively. We shall work with the constrained algorithm since we have 
the expert knowledge at oxir hands. In what follows we use Ck that goes to 0 slower 
than 0(l/fc), i.e., {l/k)/ek 0 as Jfc oo. Construct 



10 




(14) 



25 




(15) 



for estimating 9* and 




(16) 



30 



for estimating G*. 
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Averaging in both Iterates and Observations 

An alternative of the iterate averaging calls for the use of averaging in both iterates 
and outputs (observations) 



(17) 



k + l 



9k + 



k + l 



and similarly, 



k-l 



(18) 



1=0 



©A+l = ©As — 



k + l 



6* + 



1 



fc + 1 



Convergence and Rates of Convergence 

In this section, we study asymptotic properties induing convergence and rates of 
convergence of the recursive algorithms proposed in the last section. In fact, our 
attention will be devoted to the constrained algorithms using iterate averaging and/or 
averaging of both iterates and observations. Sufficient conditions will be provided, 
and existing results will be used whenever possible. The main ideas of development 
will be presented although verbatim proofs will be omitted. To proceed, we need the 
following assumptions. 

(Al) The sequences {(xjb,^jb)} and {(yib,**)} are stationary martingale difference 
processes. That is, with probability one (w.pl), 

^[i^ky <l>k)\{ock^u ^A!"i), • ' • > (^o> ^o)] = 0, and 

Elivky ^k)\iyk-u *^-i), • . . , (yo, *o)] = 0. 

(A2) The following moment conditions hold. 

EMk = J5$Jk*& = B4>kXk = B^, E^kVk = B*, 
where A4, and are symmetric and positive definite. Moreover, 

+ + 12/^1" + M < 00. 
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Remark .1. A special case of the martingale assumption deals with sequences of 
independent and identically distributed random vairables. However, more generally 
a martingale difference sequence is not necessarily independent but is uncorrelated. 
We note that allowing correlated signals is important since in the anethesia infusion 
5 problem we consider due to the operating conditionss, the influence of various noise 
sources in the enviroment, the signals generally will be correlated. The uncorrelated 
condition is only a crude approximation. We choose such a condition since it is easier 
to understand and more accessible to wider audience. It seems to be more instructive 
to present the main ideas than searching for the weakest conditions. 
10 More general correlated noise such as ^-mixing processes, which in fact allows 

diminishing correlations between the remote past and distant future; interested reader 
is referred to [5, 20] and the references therein. However, for ease of presentation, we 
confine our attention to the uncorrelated signals. To make the paper more accessible 
and appealing to wider range of audience, the assumptions given are much stronger 
15 than necessary. Our hope is that we will be able to get the main ideas acrossed 
without undue much of the technical difficulties. 

Convergence 

In what follows, we illustrate how the algorithms can be analyzed by using (17), The 
main ideas are outlined, but the verbatim proofs will be onutted. To analyze the 
20 algorithm, note that the first equation in (17) can be rewritten as 

§fc+l = ^fc + BkM^k - <l>kOk] + SkR^^k, 

where 

is the vector having shortest Euclidean length needed to bring ek(t>k[xk - (l>kdk] back 
25 ' to the constraint set Qi. 

- To prove the convergence of the recursive algorithms, in lieu of considering the 
discrete resurvie directly, we take a contiunous-time interpolation and work on ap- 
propriate function spaces. To be more specific, define 

<o = 0, and <fe = y^gj,. 

30 • 0<4<t<W 

and 

0'(*) = 5fc for t6[tfc,tfc+i), 
en*) = «'(< + tfc). 
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That is, 6^{t) is the piecewise constant interpolation of 6k and 9^{t) is a sequence 
obtained by shifting 0^(t) sequence that enables us to bring the "tail" of the sequence 
to the foreground. For a more detailed account on the interpolations, see [6, pp.90-92]. 
Similarly, define 

m(t)-l 

Rlit) = 0 for t < 0, Rlit) = ^ for t > 0, 

R^4t)^Rlit + tk)-'Rl{tk). for t>0 

and 

To proceed, define a set C{6) as follows. For 5 G flj, the interior of Cli^ C(6) 
contains only the zero element. For ^ € 9fii, the boundary of fii, C{6) is the convex 
cone generated by the outer normal at 9 at which 9 lies. Owing to the projection used, 
it is easily seen that {dk} is bounded uniformly, and as a result {^^(0} is a sequence 
of uniformJy bounded function. Moreover, it can be shown that {9^{')fR^{')} is 
equicontinuous in the extended sense ([6, p.72]). By virtue of an extended version of 
the Ascoli-Arzeld theorem, there is a subsequence that converges to some continous 
limit uniformly on each bounded time interval. Thus, we arrive the following result; 
the detailed proof may be found in [6, Chapter 5]. 

Lemma !• Under conditions (Al) and (A2), {9^{')>R!^{')} is uniformly bounded 
and equicontinuous in the extended sense. There exists a convergent subsequence with 
.limit (^(Oj-^CO) satisfying the projected ordinary differential equation 

9^B4,--A4,9 + r^, r^(t) e -C(e(t)), (19) 

where = /q Likewise, {Q^{'),R%{')} has convergent subsequence with 

limit (6(-),i2^(0) satisfying the projected ODE 

e = - A*e + r*, r*(t) e -C(e(t)), (20) 
where R^{t) = /g r*(5)d5. 

Remark 2. The significance of the (19) is that its stationary points are precisely the 
values of the parameter we are looking for. In view of the expression (19), r^(-) is 
the term needed to keep ${') in Cli. Interested reader is referred to [6, Section 4.3] for 
futher details. The significance of the r^(-) for the drug infusion process is that the 
expert knowledge allows us to keep the dynamics within a bounded region. 
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Theorem 3. In addition to the conditions of Lemma 1, suppose that there is a unique 
stationary point 6* to the ODE (19). Then dk converges to 9* w.p.l. Moreover, 
9k 9* w.pA, Similarly, suppose that there is a unique stationary point 0* to the 
ODE (20). Then Qk and 0^ converge to 0* w.p.l. 

The proof of convergence of {9^} follows from that of the argument in [6, Theorem 
5.1], whereas the convergence of is a consequence of a familiar fact in analysis 
(namely, if a sequence converges to a limit then so does its arithmetic average). 

The uniqueness assumption on the stationary point is guranteed by the choosing 
the projection bounds to be large enough such that £i < 0 and 6i > 0, i < 2n + 1. An 
argument using a form of the Kuhn-Tucker points and the ideas of active constraints 
to verify the assertion is in [6, Section 9.4]. When 9 G JlJ, = 0. Note that interior 
to Hi, the dynamics of the trajectories are determined by the ODE ^ = — A^9, 
This is a linear ODE and hence has a unique solutino for each initial condition. The 
stationary point of this ODE is precisely 0* — A^^B^. In the actual calculation 
using the recursive algorithm, if the iterates repeatedly hover the bounding surface, 
we would increase the^size of the projection region. Similar discussion applies to the 
recursive estimates 



Rate of Convergence 

For the study of rate of convergence, we concentrate on both iterate averaing al- 
gorithms and algorithms with averaging in both iterates and observations in what 
follows. 

Let us begin with (13), under the assumptions of Theorem 3, ^fc 9* w.p.l and 
^* G nj. Thus effectively, the truncation or projection can be dropped in the rate of 
convergence analysis. Define 

Vk = *fc(l/fc - *;0'). 0* = 0jb - 0*. 
Dropping the projection term and setting 



we have ^ ^ 

^fc+i =^ (/ - €kA4,)9k -f ek^k + ek[A^ - Mk% 

_ k k _ 

Using (21), we will be able to obtain a distributional limit. The statement is in 
the following theorem. Its proof can be obtained as in [6, Chapter 10). In fact, in 



(21) 
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the reference, we dealt with stochastic processes aspect of the problem by taking 
continuous-time interpolations and worked with interpolated processes, and proved a 
limit stochastic differential equation holding for the centered and rescaled sequence. 

Theorem 4 Under the conditions of Theorem 3, for the algorithm given in (13); 
(9k — d*)/V^k conver^e^ in distribution to 7V(0,E^) a normal random variable with 
mean 0 and covriance = Eii^[, and (0* — ©*)/v^jb converges in distribution to 
N{0, E^), where = Erjitf^. 

The scaling factor l/x/i^ together with the asymptotic covariances E^ and E<^ 
gives thedesired rates of convergence. We next demonstrate that iterate averaging 
and averaging in both iterates and observations lead to improved convergence speed. 
To fix the idea, we choose ek = (!/(& + 1)'^) with 1/2 < 7 < 1. Define 




Using (21), we can show that 

where o(l) — > 0 in probabilily as A; 00. Similar result holds for Wfe. We finally 
arrive at: 

Theorem 5. Under the conditions of Theorem 4, for the iterative averaging algo- 
rithms, y/kffik — 9*) converges in distribution to N{O^A'^^^^A^^) and Vki^k — ©*) 
converges in distribution to iV(0, AJ^S^AJ^). The conclusions continue to hold for 
algorisms (17) and (18). 

Discussion on Asymptotic Optimality 

Consider 0k given by (13). In view of Theorem 4, if we take ek = 0{l/k'^), with 0 < 
7 < 1. Apparenetly, the best possible scaling factor is by using Cfc = 0(l/fc). However, 
now the algorithm is stochastic. The convergence rate is determined not only by 
the scaling factor but also by its variations (variances). To examine further, take 
Bk = T/k^ where T is a matrix. According to Theorem 4, A/A(flfc-5*) is asymptotically 
normal and its asymptotic covariance in fact depends on T, i.e., = S^(r). 
Minimiring the covariance or its trace with respect to T, we obtain that the "smallest" 
covariance is given by A^^S^(>1^^)', This suggests that one may wish to construct 
another sequence to estimate E^ and use that in algorithm (13). Nevertheless, a 
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moment of reflection shows that such an idea is not feasible since it is computational 
intersive. The results in Theorem 5 show that by using averaging approaches, we 
obtain asymptotically optimal performance of the algorithms without carrying out 
the tedious estimation procedures. In addition, the use of large step size (larger than 
5 0{l/k)) provides better transient behavior. 

Simulation and Illustration 

Utility of the recursive algorithms presented in the above subsections will be illus- 
trated here by a numerical example. They will be applied to the anesthesia modeling 
in the subsequent sections. The example plant is expressed by the following Wiener 

10 system x{k) = aMk)^<^u{k-l)+w, y{k) = 2x(A:)+Ci2^g^+C22i|^g^+ 
V, where the true parameters are ai = 2; 02 = 1; Ci = 4.1; C2 = 3.5» The testing 
points for the expert information on the nonlinear function are ^1 = 0; ^1 = 5; ^3 = 10. 
For anesthesia applications, we always have y = 0 when x == 0. This information is 
already incorporated in the expression. The nominal function y = 2a; is known and 

15 contained in the expression. Prior information on unknown parameters are given by: 
range of ai = [0, 3]; range of a2 = [0, 2]; range of Ci = [2, 6]; range of C2 = [1, 5]. They 
will form fii and ^2 by Cartesian products, w and v are ii.d. uniformly distributed 
random processes, in the range [—1,5,1.5]. 

The two-step structure and the stochastic approximation algorithms are used to 

29 ■ identify 6 = [01,02]' and 0 = [C7i,C72]'. As expected, the accuracy and convergence 
of the estimates depend on the input signal (for richness in its probing capability), 
step sizes, and prior information. Many variations on these conditions are tested: 
for inputs: the step input, periodic inputs of different periods and typical anesthesia 
infusion profiles; for step sizes: constant sizes of different values, variable sizes of 

25 diflFerent initial values and decaying rates; for prior information: different ranges. A 
typical estimation error trajectory is illustrated in Figure 3, which has a periodic 
input and variable step size efc = 0J5ft"'°*^ for both 9 and 9 estimators. The resulting 
estimation: 9 = [2.04, 1.03]' and 6 = [4.21, 2.9951]', much improved accuracy from 
the prior expert assessment. 

30 Model Development for Patient BIS Responses 
to Drug Infusion 

It is a great challenge to characterize mathematically a patient's response to drug in- 
fusion. Due to large deviations in physical conditions, ages, metabolism, pre-existing 
medical conditions, and surgical procedures, patients demonstrate high nonUnearity 
35 and a large variations in their responses to drug infusion. Also, detailed models of 
clinical subjects on their physilogical processes are inherently very complex and not 
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suitable for clinical applications. On the other hand, anesthesiologists have been 
administering drug infusion successfully with only limited information on patients, 
such as weights and medical conditions. Control strategies of an experienced anes- 
thesiologist are based intuitively on basic characteristics, such as how sensitive the 
5 patient is to a drug infusion. The modeling approach, previously introduced captiures 
this imderstanding and is briefly summarized below. 

Model Structures 

The model contains three building blocks, as shown in Fig . 1 The first two blocks 
represent the patient response to popofol titration and bolus injection, respectively. 
10 The tlurd block characterizes the patient response to predictable stimulation such as 
incision, closing, etc. For concreteness, we will assume that propofol is the anesthesia 
drug although the model structure is generic and valid for other anesthesia drugs. 
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Propofol titration is administered by an infusion pump, whose dynamics is rep- 
resented by a first-order dynamic delay where Ti is the time constant of the 
infusion pimap. Ti typically ranges from 0.2 to 2 seconds and is determined from 
experiment data. 

5 The key component of the model is the patient dynamics, which is defined by three 
basic elements: (1) Initial time delay Tp after drug infusion, represented by a transfer 
function e~'>*; (2) Time constant Tp representing the speed of response, ^r^; (3) 
A nonlinear static function representing sensitivity of the patient to a drug dosage 
at steady state, represented by a nonlinear fxmction /p. The meanings of these three 

10 components are depicted in Figs. 2a and 2b. 

Two possible model structures tor titration are shown in Fig. 3. depending on 
if the nonlinear function preceeds the linear dynamic system (Hammerstein Model) 
or follows it (Weiner Model). Both model structures are viable in this application. 
For simplicity, we shall focus on the Wiener model. The delay time Tp, time constant 

IS Tpj and the sensitivity function vary greatly among patients and must be established 
individually for each patient. 

Bolus injection inputs a large amount of drug in a short time interval. The bo- 
lus model has a structure similar to the titration model, but with different model 
parameters, that are correlated to those of the titration model. 
20 One of the most important control tasks in infusion control is to predict the im- 

pact of surgical stimulation on the BIS levels. Surgical stimulation is the main cause 
that induces large fluctuations in anesthetic df^^^th during a surgical procedure. Early 
prediction of the impact will allow advanced drug adjustment before a feedback 
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signal from the BIS monitor activates control actions. Typical surgical stimulation 
includes initial surgery preparation, incision, closing, etc. Significant surgical stimu- 
lation is predictable in advance and can be roughly characterized by their severity in 
affecting BIS values* Here we use classification levels from 1 to 5 to represent severity 
of stimulation: 1 being a minor stimulation and 5 the most significant stimulation. 
Patient responses to surgical stimulation are represented by the third model block. 
The block contains also a time delay and a transfer function to reflect its 
dynamics and a lookup table for its impact on the BIS levels. Since a patient's BIS 
sensitivity to stimulation depends on hypnosis state, a scaling factor a is introduced. 

Stimulation Level x a = Scaled Stimulation Level 

The scaled stimulation is then used to predict the impact on the BIS measurements. 
Information on an upcoming stimulation will be used as a pre-warning to provide a 
predictive control so that the impact of the stimulation can be promptly compensated 
by early tuning of drug infusion. 

The BIS monitor measures the BEG signals and performs bi-spectrum statistical 
analysis. Due to the processing time and data windows required in this procedure, 
the BIS monitor introduces a pure time delay (from processing time) and a dynamics 
delay (related to the size of data windows). Its relevant dynamics to control can be 
represented by a system transfer function e^''^^ cascaded with a nonlinear limiter 
with lower limit 0 and upper limit 100 (defining the BIS range). 
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There are many sources of disturbances that will also affect the BIS readings. 
These include unpredictable stimulation to the patient, other drugs, electrical noises, 
sensor attachment movements, environment noises, etc. These disturbances are un- 
predictable and cannot be compensated by early predictive controL Their impact on 
BIS levels will be minimized by a feedback action. 

Model Verification and Development by Clinical Data 

To verify the utility of the model structure and develop individual patient models, 
clinical data are collected, after an approval from the IRB board. One of these data 
sets is described below to illustrate the process and results. 

The anesthesia process was administered manually by an anesthesiologist and 
lasted about 76 minutes, starting from the initial drug administration and continiung 
until the last dose of administration. The trajectories of titration and bolus injec- 
tion, as well as the corresponding BIS values are shown in .Fig. 5. Scaled surgical 
stimulation trajectories are determined by the anesthesiologist and shown in Figure 
6. 

The data were then used to determine model parameters and function forms. The 
Wiener model is used in this case. The actual BIS response is compared to the model 
response over the entire surgical procedure. Comparison results are shown in Figure 
30.Here, the inputs of titration and bolus are the recorded real-time data. The model 
output represents the patient response very well. In particular, the model captures 
the key trends and magnitudes of the BIS variations in the surgical procedure. This 
indicates that the model structure contains sufficient freedom in representing the 
main features of the patient response. Also, the impact of surgical stimulation is 
captured by the model. These features are of essential importance in designing control 
strategies for computer-aided drug infusion. 

Real-Time Identification of the Patient Model 

In the previous section, by using an off-line modeling process we have shown that 
the model structures " can represent satisfactorily the patient 

response to drug infusion, in cumcal applications, this model must be developed 
individually for each patient and modified over the surgical procedure. Consequently, 
fast real-time modeling becomes inevitable. The main goal of real-time modeling is 
to obtain the model parameters quickly and accurately. 
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Employing the identification method and the patient model 

we will develop a fast identification method for patient models in anesthesia 
infusion. The algorithm contains two modules: Expert Classification Module (ECM) 
and Real-Time Learning Module (RTLM). For simplicity, we shall focus on the titra- 
5 tion model. The Wiener model entails a linear time-invariant system (time delay and 
time constant) followed by a memoryless nonlinear function (sensitivity function). 
Mathematically, this relationship is expressed as 

where U{3) is the Laplace transform of the input u (the propofol titration rate). The 
10 BIS value is related to y by 

BISjb = 100 - Sat(yA) 

where Sat is the saturation function defined by S^t{y) = y, if 0 < y < 100; Sat(y) = 0, 
if y < 0; and Sat(y) = 100, if y > 100, This expression simply represents the fact 
15 that the drug infusion will reduce the BIS value from its maximum value 100 (fully 
awake). 

Since the model parameters Tp and Tp have clear physical meaning and are com- 
patible with the accessment of an anesthesiologist, this model allows direct inputs 
from expert knowledge. For instance, a patient with fast drug response will have a 
20 smaller Tp and Tp. It should be emphasized that such a subjective accessment is not 
accurate and can serve only as initial parameter ranges. 

The sensitivity function fp is a monotone increasing function (the higher the drug 
infusion rate, the stronger the drug impact), and fp{0) = 0. In this application, m = 3 
is sufiicient and typical set points are selected to be xi = 50, X2 = 100, = 150. To 
25 accommodate the condition fp{0) = 0, the base functions in (1) is slightly modified 
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X\ X2 '*'3 

The coefficients difer significantly from patient to patient. Determination of the 
coefficients Ci, C2, C3, and Tp and Tp, will be the main task for real-time learning- 

5 The Expert Classification Module 

Since control strategies must be decided at the outset of the surgery, one must obtain 
initial patient characteristics to determine the initial control action. Without real- 
time data, these parameters must be initiated by expert assessment. The expert 
knowledge of an anesthesiologist on a patient's drug response pattern is obtained 

10 from a patient*s medical history and conditions. An anesthesiologist employs the 
inforraatioin to give an initial approximate accessment, usually expressed in linguistic 
terms such as "this is a patient who may be very sensitive to propofol." 

The Expert Classification Module (ECM) represents a mapping from the linguistic 
assessment to initial values Tp, Tp, and Ci,C72,C3. Here, we ask an anesthesiologist 

15 to assign a relative value between 1 and 10 to represent a patient's resp[onse to drug 
infusion, where 1 represents the slowest and 10 the fastest. A computerized mapping 
translates the assessment to a bounded parameter set on Tp and Tpy as shown in Figure 
10. 

For an individual patient, the anesthesiologist provides an initial assessment on 
20 drug sensitivity: Given the titration rate Xi, the corresponding estimated raiiges of 
BIS values [y.,yi], < = lj2,3. This leads to the prior parameter uncertainty set 
^^2 = (EijI/i] ® [y^^y2] ® [Us^Vzh where ® is the Cartisian product. 

This framework allows seamless incorporation of expert knowledge into the math- 
ematical model structure even when this knowledge is partial and approximate. 

25 The Real-Time Learning Module 

Starting from the initial model uncertainty sets Qi and Shi ^be Real-Time Learning 
Module (RTLM) employs real-time data to learn the patient characteristics ajid to 
improve model accuracy. Algorithm 1 is employed in this module to identify the 
parameters recursively. 

30 Verification 

The algorithm is evaluated by using the collected clinical data. Under manual control 
of drug infusion by an anesthesiologist, the popofol titration, bolus injection, surgical 
stimulation, and the patient response are recorded. The actual samling interval for 
the BIS trajectory is 5 seconds and the rate change of the drug infusion is timed 
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and recorded manually. For the convenience of data processing, the data are then 
reformated and resampled (interpolation) to form a synchronized input-output data 
with sampling interval 1 second. 

To illustrate, we shall concentrate on identification of the titration model. To 
5 estimate the parameters in the titration block, the data in the interval where the 
bolus and stimulation impact is minimal are used. A data window of 8 minutes is 
used in the identification process. The identification proceeds as follows. 

First, a grid of potential values of the time delay Td and the time constant Tp is 
generated. In this particular application, a grid width of 2 seconds is used for Td and 
10 a grid width of 5 seconds is used for Tp, resulting in a grid space G of size x nip, 
(Tj, r/), i = 1, - . . , rrid, j = 1, . . . , mp. Due to robustness in feedback and predictive 
control, this level of precision is sufiBLdent. For each selection of the values from the 
grids, we estimate the sensitivity function l^ optimized data fitting (least-squares). 
Then the fitting errors are compared over the grid space G. The optimal (S^^STp^*) 
is selected that minimizes the fitting error. 

The actual BIS response is then compared to the model response over the entire 
surgical procedure. Figure 31 shows the model BIS trajectory versus the actual BIS 
data from the patient. 
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Real-time medical decisions are exemplified by general anesthesia for attaining 
an adequate anesthetic depth (consciousness level of a patient), analgesia (pain 
management), sedation control in ICU (intensive care units), fluid resuscitation in 
trauma cases, circulation control (ventilation and mechanical circulation), etc. In all 
these medical decision processes, one of the most critical requirements is to predict the 
impact of the inputs (drag infusion rates, fluid flow rates, airway pressure and flow rates, 
etc^ on the outcomes (consciousness levels, pain scores, blood pressures, heart rates, 
oxygen saturation, efc.)- This prediction capability is not necessarily "control-oriented,"* 
although its control implication is obvious. It can be used for display, warning, 
predictive diagnosis, decision analysis, outcome comparison, etc. 

The core function of this prediction capability is embedded in establishing, in 
real-time, a reliable model that relates the (multiple) inputs to the (multiple) outcomes. 
Apparently^ this is a real-time identification problem. This problem offers a great 
challenge and opportunity for advancing system identification. 

The major challenges that must be addressed include: 

(1) Reliability and safety are mandatory. Theoretical analysis and limited 
simulation must be further enhanced by a more comprehensive evaluation process. 
Eventually reliability must be established statistically via clinical trials. 

(2) The characteristics of patient responses demonstrate significant nonlinearity 
and time variation. 

(3) Unlike industry applications, patient responses depend critically on patient 
medical conditions, surgical procedures, and drag interactions; and hence they are not 
repeatable. Models must be established individually and in real-time. 

(4) This is not a data rich environment. At the starting time interval before 
substantial data become available one must rely on expert assessment to initiate the 
model and decision process. This implies that the modeling process must allow a 
seamless combination of assessment and data-generated estimation. 

(5) Since data points are limited, low-complexity models and fast identification 
are always preferred. 

(6) Due to medical constraints, the inputs are not allowed to be arbitrarily 
selected for best identification experiments or probing excitation. 

These challenges, together with traditional focus of system identification on 
linear, well defined, and repeatable environment, have resulted in a gap between 
theoretical results in identification and applications to medical applications. 
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Patient Model Structures 

In an attempt to understand the potential of system identification in diis 
important application area, we have identified anesthesia decisions as a typical setting 
for studying the system identification problem. Some preliminary findings on the utility 
of Wiener model structures (Hammerstein models can be similarly applied), 
combination of assessment and real-time estimation, and identification algorithms are 
summarized below. 

The goals of general anesthesia are to achieve hypnosis and analgesia 
simultaneously throughout a surgical operation while maintaining the vital functions of 
the body. One of the most critical tasks of anesthesia is to attain an adequate anesthetic 
depth. 

The newly marketed BIS monitor provides a viable direct measurement of 
anesthesia depth. To facilitate control strategy development for improved anesthesia 
performance, it is highly desirable to develop representative models of BIS responses 
to dmg infusion. Since patients differ dramatically in metabolism, pre-existing medical 
conditions, and surgical procedures, individualized models must be established for each 
patient with limited patient information and data points. A knowledge-based and 
control-oriented modeling approach was introduced recently and applied to develop a 
feedback and predictive control strategy for anesthesia infusion. The models retain only 
the key characteristics of patient responses that are essential for control strategy 
development and can be determined by either data or expert knowledge. Furthermore, 
an identification algorithm has been introduced for updating the patient model in real- 
time. 

Due to complications in human metabolism and nerve systems, we do not 
perceive that the system dictates any specific model structure as especially superior. 
Simple model structures that can capture the key characteristics and maintain 
physiological meanings of model parameters will be preferable since they will allow 
easy interface with physicians and fast model updates. 

As previously noted, we have tested the validity of Wiener structures in 
representing the patient dynamics. The patient dynamics is defined by three basic 
elements: 

(1) Initial time delay after drug infusion; 

(2) Time constant representing the speed of response; 

(3) A nonlinear static function)^ representing sensitivity of the patient to a dmg 
dosage at steady state. 
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The meanings of these thiee components are depicted in Figs. 2a and 2b. 

This model structure can be mathematically represented by a Wiener model with 
input u and output y, x=Gu H-^w, y=f(x)+"v where G is a dynamic linear system, /is a 
memoryless nonlinear function, and "w and "v are noises, is an intermediate variable 
that cannot be directly measured. In our applications here, G and / represent the 
transient and steady-state behavior of a patient's response to a drug infusion input, 
respectively. G can be discretized and represented by an ARMA model: 

where 

or approTumated by a MA model with sufiicientiy high orders. 

The sensitivity function/can be modeled by either a linear combination of some 
pre-selected base functions or a parameterized nonlinear function. The delay time ^ , 
time constant T^, and the sensitivity function/ vary greatiy among patients and must be 
established individually for each patient. 

To verify the utility of the model structure and develop individual patient 
models, clinical data are collected. One of these data sets is described below to illustrate 
the process and results. The anesthesia process was administered manually by an 
anesthesiologist and lasted about 76 minutes, starting from the initial dmg 
administration and continuing until the last dose of administration. 

To verify the capability of the model stracture, the data were used to determine 
model parameters and function forms, by an off-line estimation method (an optimal 
nonlinear LS estimator). The actual response is compared with the model response over 
the entire surgical procedure. Comparison results are shown in Fig. 10a to lOd. Here, 
the inputs of drug titration (continuous drug flow) and bolus (drug injection of a short 
duration) are the recorded real-time data. The model captures the key trends and 
magnitudes of the output variations in the surgical procedure. This indicates that the 
model structure, though very simple, contains sufficient freedom in representing the 
main features of the patient response. Also, the impact of surgical stimulation is 
captured by the model. 
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Reliability Evaluation of Identification Algorithms 

Despite well established theoretical results on many identification algorithms, our 
studies indicate that to enhance reliability of system identification in medical 
applications, some well developed algorithms must be significantly re-evaluated and 
modified beyond routine convergence analysis. Our criteria for evaluating an algorithm 
include: 

(1) Accuracy and convergence (estimation errors); 

(2) Robustness (error excursion frequencies); 

(3) Model complexity (number of parameters); 

(4) Time complexity (data points required); and 

(5) Usage convenience (requirements of delicate tuning or sensitive dependence 
on initial values are not desirable). 

To understand the reliability of various identification methods on Wiener models 
and anesthesia patient models, the following platforms are used as a benchmark. 
Platform 1: A Simulated Wiener System 
The example plant is expressed by the Wiener system: 



where the true parameters are a/=2; a2=U Cy=4.1; C^3,5. The testing points for the 
expert information on the nonlinear function are ^=5; ^=10. For anesthesia 
appUcations, we always have ypO when ^^^=0. The nominal function yp2x, is known and 
contained in the expression. Prior information on unknown parameters are given by: 
range of a;=[0, 3]; range of a2=[0, 2]; range of Cy=[2, 6]; range of C2=[l, 5]. 
w, and V, are i.i.d. uniformly distributed random processes, in the range [-1.5, 1.5]. 
Platform 2: An Anesthesia Patient Model 

From the patient data illustrated in Figure \ref { model } , the titration model is first 
extracted by an off-line method. The method eliminates the impact of surgical 
stimulation and drug impact fi:om bolus injection, and produces a BIS response of the 
dmg propofol titration on the patient as shown in Figs. 19a and 19b. 
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Recursive Algorithms 

Model parameter updating algorithms must be used to estimate the parameters 
of the Wiener model, especially the linear dynamic part. There have been established 
many valid recursive algorithms in a wide variety of applications. Our objective here 
is to evaluate the accuracy and convergence of different recursive algorithms under this 
application* Here a linear system is used 



where c is an i.i.d. Gaussian disturbance with zero mean and unit variance. 

The evaluated algorithms include the following four types. For the first three 
algorithms, since the tunable parameters have dramatic impact on their performance, 
their values are optimized in each run. 
(1) Adaptive FUtering 




(2) Adaptive Filtering with Averaging 



(Z>,^4),0<r<l; 



1 




1 



t+l 



(3) One-step Modified Optimal Projection 
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(4) Recursive Least Squares 




The algorithms are started with the same initial estimate, which is obtained by 
using an LS estimate with 10 data points. Then the recursive algorithms are perforaied 
for 1000 data points. For each algorithm, the estimation errors during the last 20 time 
points are averaged. The averaged error is used as a measure of identification accuracy. 
The simulation is repeated 50 times. The 50 estimation errors are plotted in the right 
column of Figs 20a to 20h for algorithm comparison. Also, the 50th error trajectory for 
each algorithm is illustrated in the left column. 

It seems clear that the recursive least squares algorithm outperforms the rest in 
two aspects: (1) It provides more accurate estimates; (2) It does not require any variable 
tuning. It is noted that we have optimized variables for each run in other algorithms. 
In real applications, these variables cannot be optimized over each sample path. 
Consequently, the actual performance of these algorithms will be worse than what we 
have demonstrated in this simulation. In fact, when the variables are fixed, we 
encountered occasionally parameter divergence in adaptive filtering or projection 
algoritiims. In general, the RLS may require slightly more computational time. 
However, this is not an issue in this case since the model order is small in patient 
models. Consequently, we will use LS algorithms in our identification of patient 
models. 

Block Recursion 

A two-step recursive estimation algorithm has been proposed for identifying the patient 
model in real-time. The procedure can be briefly summarized as follows; 
Two-step Procedure. We start at time fc=0 with 



(the parameters for the linear part) and 



@Q (the parameters for the nonlinear part). 
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Step 1: At time fc, the available 0 q is used to map the output y^^ backwards to obtain 

/\ 

the "intemiediate" output . Then and Jc^ are used to update ^^^j . 
Step 2: From the obtained j we map u^^+y forward to obtain x^^^^ , Then Xj^^^ 

and yjt^; are used to update from . 

5 Selection of recursive algorithms, averaging, projection into a bounded set, step 

size selection, and convergence analysis have previously been presented. This algorithm 
is computationally very efficient. 

To evaluate its reliability, a known system is used that works reasonably well 
when initial estimates are close to true values. However, when initial estimates are 

10 furth^ away from true values the algorithm demonstrated irregular behavior that 

includes parameter div^gence, slow parameter drifting, and large output prediction 
errors. In other words, it is not a reliable algorithm. 

Consequently, the algorithm is modified into a block recursive algorithm. In this 
modified algorithm, the parameters are not updated upon receiving each observation. 

15 But rather, a block of m observations is used collectively to update the model 

parameters. Starting with the most recent estimate as an initial condition, the algorithm 
searches for the optimal parameters that minimize the output errors between the model 
output and measured output in the data block. This modified algorithm have shown a 
greatly improved reliability. For the system (2), we used the block size 10 to. run a 

20 simulation of 800 data points with a random input. The simulation was repeated 50 

times and demonstrated good prediction capability in each ran. A typical result is shown 
in Fig. 21. On the other hand, parameter convergence is not always demonstrated. 
Optimization J^fficiency 

For a nonlinear system, optimization carries potentially large computational 

25 burdens. As shown in the previous section, overly simplified algorithms may suffer loss 

of reliability. In search of efficient optimization methods during a block recursion, an 
embedded optimization method was employed. In this method, for given parameters of 
the nonlinear part of the model, the parameters of the linear part are optimized by the 
least squares algorithm that is computationally very simple. Consequently, only the 

30 parameters of the nonlinear function must be searched. In our evaluation, global grid 

search, local grid search (around the initial estimate), simplex optimization methods 
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were tested. The simplex method finds similar parameter values but requires on average 
only about one twentieth of computational time required by the grid search methods. 
Parameterization of Nonlinear Functions 

The memoryless function / in the Wiener model represents steady-state 
sensitivity: The impact of the input (dmg amount) on the outcome. Theoretically, a 
memoryless and continuous nonlinear function can be approximated by a polynomial 
function to any degree of accuracy when the order of the polynomial is allowed to 
increase. However, the parameters of the polynomial usually do not carry any physical 
meaning. It has been proposed that the following functions are more suitable to connect 
with expert assessment. 

A physician's knowledge can be expressed as follows: Given m typical x 
amounts ^i, ^n,, the approximate ranges of the predicted steady-state outcomes are 

iy^^yii^-^'^ly^^ym^ • Weparameterlze/by: 
where ®i(;c) is the m* order polynomial: 



<&i(jc) satisfies *i(5j)=0, j ^ i; Consequently, the prior expert information 

becomes precisely the prior interval information on the parameters 0=[Ci,...,CJ E 

[y.yyi] ,i=^ 

This knowledge will serve as the prior information on the unknown parameters 
0 and ©• Consequently, the parameters will have clear physical meanings and expert 
assessment can be integrated seamlessly into this parameterization. 

Interestingly, for computational purposes this parameterization implies further 
complications. Understanding that the nonlinear function in this application represents 
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a drag's steady-state impact on the patient outcome. Hence it must obey the basic 
monotone principle: the more the drag infusion rate, the more the impact on the patient. 
Unfortunately, this requirement is not naturally * embedded in the polynomial 
parameterization. In fact, in a box of parameter values for the polynomials in the 
5 resulting nonlinear functions are not always monotone. Also, for the same parameters, 

the corresponding function may become non-monotone when the range of the drag rates 
is enlarged. Consequently, in search of optimal parameters during system identification, 
monotonicity must be laboriously checked in each evaluation and sonjietimes leads to the 
failure of a local search. 

10 A remedy of this situation is found by using a different parameterization of the 

nonlinear function. The following two-parameter function is employed: Suppose that 
X takes values in [0,Z>]. 



where erf (.) is the standard error function, and the parameters r and a define variation 

15 in y at jcs=fc and curvatures of the erf function in jc e [o,fc] . Since it is very simple to find 

a box in r-^space under which the resulting nonlinear functions are always monotone, 
it greatly improves time complexity in search algorithms. In Figs. 24, 25, and 26, the 
ranges of searched nonlinear functions are included in the bottom plots. 
Model Complexity and Time Complexity 

20 Model complexity and time complexity are especially important in medical 

applications due to limited data points, input probing richness, requirements of prompt 
decisions, and time varying natures of patient models. In this evaluation the patient 
model in Figs. 19a and 19b is used. Since actual model order and stractures are 
unknown to the designer, we simply use either MA or ARMA models of various orders 

25 and evaluate their prediction capability. The input data are the real drag infusion rates 

which do not provide much richness in excitation. The total surgery lasted about 4000 
seconds. We used the first 400 seconds of the data for system identification for ARMA 
models, and 1000 data points for MA models since their orders are much higher. Then 
the prediction capability of the identified model is evaluated by applying it to the entire 

30 data. It is noted that during the first 400 seconds, the input changes its value only once. 

Conceptually, it is arguable that for a linear stable system, it is always possible 
to approximate it by a FIR model (moving average). This model stracture can 
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significantly simplify theoretical analysis. However, a comparison of Figs. 22 and 23 
shows clearly that to capture the dynamics of the actual system, the order of the MA 
model must be very high. This is highly undesirable in this application since 200 
parameters in this model is a very high model complexity. It also increases 
identification time complexity. Since an overall patient model contains many model 
blocks of similar types, the MA structure will imply a very large parameter set for 
system identification when the overall system is considered. 

Figs. 24, 25, and 26 show the prediction capability of ARMA models of orders 
1, 4 and 10, respectively. It is seen that ARMA models can capture the dynamics of the 
system with much lower orders than MA models. 

In sum, it is noted that: (1) Wiener or Hammerstein models have great potential 
utility in modehng patient dynamics; (2) Reliable identification is a stringent 
requirement that mandates a modified evaluation criterion on identification algorithms; 
(3) Within the fi:amework of reliable identification, all components of a system 
identification must be carefully examined in terms of robustness, reliability, model 
complexity, time complexity, accuracy, and tracking capability for time varying 
parameters. 

Although the invention has been described in terms of specific embodiments and 
applications, persons skilled in the art may, in light of this teaching, generate additional 
embodiments without exceeding the scope or departing from the spirit of the claimed 
invention. Accordingly, it is to be understood that the drawing and description in this 
disclosure are proffered to facilitate comprehension of the invention and should not be 
construed to limit the scope thereof. 
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